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ABSTRACT 

The ionising continuum from active galactic nuclei (AGN) is fundamental for inter¬ 
preting their broad emission lines and understanding their impact on the surrounding 
gas. Furthermore, it provides hints on how matter accretes onto supermassive black 
holes. Using HST s Wide Field Camera 3 we have constructed the first stacked ultravi¬ 
olet (rest-frame wavelengths 600 - 2500 A ) spectrum of 53 luminous quasars at 2 ~ 2 . 4 , 
with a state-of-the-art correction for the intervening Lyman forest and Lyman con¬ 
tinuum absorption. The continuum slope (/„ cx !/“■') of the full sample shows a break 
at ^912 A with spectral index a„ = — 0.61 ± 0.01 at A > 912 A and a softening at 
shorter wavelengths (a„ = — 1.70 ± 0.61 at A ^ 912 A). Our analysis proves that a 
proper intergalactic medium absorption correction is required to establish the intrin¬ 
sic continuum emission of quasars. We interpret our average ultraviolet spectrum in 
the context of photoionisation, accretion disk models, and quasar contribution to the 
ultraviolet background. We find that observed broad line ratios are consistent with 
those predicted assuming an ionising slope of a; on =—2.0, similar to the observed ion¬ 
ising spectrum in the same wavelength range. The continuum break and softening are 
consistent with accretion disk plus X- ray corona models when black hole spin is taken 
into account. Our spectral energy distribution yields a 30 % increase to previous esti¬ 
mates of the specific quasar emissivity, such that quasars may contribute significantly 
to the total specific Lyman limit emissivity estimated from the Lya forest at z < 3 . 2 . 

Key words: accretion, accretion discs - galaxies: active - quasars: general 


1 INTRODUCTION 

Considerable effort has been devoted to characterising the 
shape of the quasar (QSO) ionising continuum in the ultra¬ 
violet (UV) over the past years. The spectral energy dis¬ 
tribution (SED) of active galactic nuclei (AGN) shows a 
prominent bump, the so-called “Big Blue Bump” (BBB), 
which appears to peak in the UV and decline at higher 


energies (Sanders et al. 1989 Elvis et al. 1994). However, 


the intrinsic position and possible luminosity dependence of 
the BBB has not been properly estimated due to the lack 
of UV observations corrected for the intergalactic medium 
(IGM) absorption by neutral hydrogen along the line of sight 
(|Richards et al.||2006a| |Trammell et al.||2007| |Shang et al. | 
2011||Elvis et al |2012| >. 

One of the main predictions of accretion disc models is 
that the thermal disc should peak at bluer wavelengths for 


smaller black hole masses (i.e., T max oc (L/I/Edd) 1 ^ 4 A7 B 1 ^ 4- 


Shakura fe Sunyaev|1973|. However, both the presence and 


the position of the break are highly dependent on the IGM 
correction considered. 

From a theoretical perspective, the AGN ionising 
continuum is crucial for interpreting broad and narrow 
emission-lines observed in AGN spectra, and their relative 


ratios (see Stern et al. 2014 and Baskin et al. 2014 for re 


cent models of these emission line regions). There is a general 
consensus that AGN emission lines are produced by a pho- 
toionising continuum (extending from optical-UV to X-ray) 
which emerges from an accretion disc around the black hole 
and by a hot (T ~ 10 8-9 K) plasma of relativistic electrons 
that Compton up-scatter the photons coming from the disc 
( Haardt fe Maraschi||199l| 19931. 

Quasars are relevant (maybe even the dominant) 
sources of ionising photons that determine the ionisation 
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state and the temperature of the 2 ~ 3 IGM (e.g. 

Haardt & 

Madau|1996, 2012; Meiksin & White][2003; Faucher-Giguere 

et al. 112009 

). While not numerous enough at 2 > 

6 to have 


contributed significantly to Hi reionisation (e.g. Meiksin 


The differences among various surveys may arise from 
several factors, such as the small number of observations 
covering A < 1216A (e.g. less than 20 AGN at 2 > 2 in 
the T02 sample), and the crucial placement of the “con- 


20051 Jiang et al. 2008; Shankar & Mathur 2007; Willott 

tinuum windows” in order to construct the quasar ionising 

et al. 2010 

Fontanot et al. ||2012 

2014 

1 , they are likely the 

continua. Although the latter may not bias significantly the 


only sources responsible for the reionisation of He 11 at 2 ~ 3 
(|Miralda-Escude et al. 2000 Faucher-Giguere et al. |2008 


Furlanetto 2009| McQuinn et al.||2009[ |Haardt fe Madau| 
2012||Compostella et al.|2013l. All these studies rely on pa- 


rameterizations of the quasar SED in the BBB region mo¬ 
tivated by existing observations, the uncertainties of which 
are rarely discussed. 

From an observational point of view, composite spec¬ 
tra of AGN were constructed by taking advantage of sev¬ 
eral surveys (LBQS,|Francis et al.|199l FIRST, Brotherton 


et al. 2001 SDSS, Vanden Berk et al. 20011. In all these 


studies, the rest-frame optical composites indicate that the 
continuum can be described by a power law of the form 
/„ oc u au , with —0.5 < a v < —0.3 in the wavelength range 
1200-4000 A. The first composite that suggests a softening 
in the far-ultraviolet (blueward of Lya) was reported by 


Zheng et al. (1997 Z97 hereafter) who analysed 101 quasars 


from the Hubble Space Telescope ( HST ) in the redshift range 
0.33 < 2 < 3.6, covering the wavelengths between 350 and 
3000 A. This softening was interpreted as Comptonization of 


the thermal disc emission in a hot corona above the disc ( Cz- 
|erny fc Elvis|[l987| , due to the similarity between the slope 
of —1.8 found by Z97 at A < 1216A and the one measured by 


Laor et al. (1997) for a sample of radio quiet X-ray selected 


quasars (a v ~ —1-7), for which simple accretion disc+X- 
ray corona models were utilised. The break would thus in¬ 
dicate the peak of the BBB. The work by Z97 has been 
extended by |Telfer et al.| ( |2002| hereafter T02) with more 
than 80 quasars from HST over a similar redshift range. 
This analysis confirmed the findings of Z97 that the UV 
spectral continuum can be parametrised by a broken power 
law with the break in the vicinity of the Lya, with a con¬ 
tinuum slope of a„ = —1.57 for the radio quiet T02 sample. 
These results are at variance with those presented by |Scott| 
et al. (2004] S04 hereafter), who considered more than 100 


AGN at 2 < 0.1 observed with the Far Ultraviolet Spectro¬ 
scopic Explorer ( FUSE ), covering the rest-frame wavelength 
range 630—llOOA. The spectral slope of the FUSE compos¬ 
ite spectrum is a„ = —0.56, significantly harder than pre¬ 
vious estimates in the far-infrared from HST studies. IShullI 
et al. 2012J S12 hereafter) have measured the AGN ion¬ 


ising continua in 22 AGN at 0.026 < 2 < 1.44 using the 
Cosmic Origins Spectrograph (COS) on HST , covering the 
rest-frame wavelength range 500— 2000A. The COS compos¬ 
ite shows a break at A ~ 1000A, in line with the previous 
estimates by Z97 and T02, but with slightly harder spectral 
index (a„ = —0.68 at A = 1200 — 2000A and a„ = —1.41 at 
A = 500 — 1000A). Recently, Stevans et al. (2014, S14 here¬ 
after) considerably improved the S12 analysis by adding 137 
AGN to the original sample for a total of 159 sources, se¬ 
lected from the COS archive at redshifts 0.001 < 2 < 1.476. 
This new COS composite is fully consistent with the one in 
S12 and shows similar spectral indexes. We note that pre¬ 
vious estimates of the break were performed with very few 
spectra (~10 in T02, 3-12 in S12, ~20 in S14) contributing 
at short wavelengths (e.g. < 700A). 


results, a possible uncertainty arises from the selection cri¬ 
teria involved in defining the various samples. The tradi¬ 
tional strategy is to extract any source available from the 
HST/FUSE archives and apply several cuts such as red¬ 
shift, signal-to-noise, and wavelength coverage. However, the 
brightest UV sources were usually targeted for ultraviolet 
spectroscopy, primarily for absorption line studies. Further¬ 
more, any archival research has the tendency to include more 
peculiar objects that were selected for special investigations, 
and thereby the same objects were re-observed with every 
generation of UV spectrographs. These issues likely bias any 
UV archival sample to be very bright in the ultraviolet with 
a highly heterogeneous selection function. 


An additional source of variance comes from the correc¬ 
tion for the IGM absorption employed by different authors. 
Although the technique adopted by T02 and S04 was simi¬ 
lar, i.e. a statistical correction for the unidentified absorbers 
by considering an empirical parametrisation (see |Petitjean| 
et al.| fl993 Dave & Tripp 20011, the overall correction ap¬ 
plied was very different as pointed out by S04. In fact, S04 
noted that the correction considered by T02 was less than 
1% over the whole rest-frame wavelength range, which un¬ 
derestimated the number of Lya absorbers by a factor of 
~ 50 at 2 = 0.1. Another simplification usually adopted in 
previous studies is the use of a single correction for quasar 
samples spanning a large range of redshift, which, because 
of strong evolution in the IGM absorption, likely results in 
significant errors. 


In this paper, we present the first average quasar spec¬ 
trum in the rest-frame wavelength range 600-2500 A cor¬ 
rected for intervening H 1 Lyman forest and continuum ab¬ 
sorption with state-of the-art IGM transmission functions 


calibrated to the most recent observations (Prochaska et al. 


2014). The sample consists of 53 2 em — 2.4 quasars from the 


HST survey for Lyman limit absorption systems (LLSs) us¬ 


ing the Wide Field Camera 3 (WFC3) presented in O’Meara 


et al. (2011 011 hereafter). The structure of this paper is 


as follows. In Section [2] we discuss the sample and the se¬ 
lection criteria. In Section [3] we describe the technique to 
construct the stacked spectrum, whilst the IGM transmis¬ 
sion curves adopted to correct the observed average spec¬ 
trum are presented in Section [4] where we also describe our 
IGM corrected stack with uncertainties. The results and im¬ 
plications of our analysis are discussed in Section [5] while 
the conclusions are presented in Section [6] 


We adopt a concordance flat A-cosmolo gy wit h Hg = 
70kins -1 Mpc -1 , = 0.3, and Ha = 0.7 (Komatsu et al. 


2009). Unless noted otherwise, we will distinguish between 


the following wavelength ranges in the UV: (i) the near UV 
(NUV; 2000-3000 A), (ii) the far UV (FUV; 912-2000 A), 
and (iii) the extreme UV (EUV; A < 912 A). 
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2 THE DATA SET 


The data sample employed in the present analysis comes 
from a survey performed with HST using the WFC3 in¬ 
strument. Quasar sample, HST observations, and reduction 
procedures are described in detail in Oil. In this section we 
provide a brief summary of this data-set. 

The Oil survey consists of 53 quasars selected from 
SDSS Data Release 5 ( Schneider et al.|2007 ( with g < 18.5 
mag, 2.3 < z em < 2.6, observed with the WFC3/UVIS- 
G280 grism in Cycle 17. These data were taken specifi¬ 
cally for the scientific goal of surveying the abundance of 
strong Hi Lyman limit absorption features at 1.2 < z < 
2.5. Flux and wavelength calibrated ID spectra for each 
quasar in this sample are extracted using customized soft¬ 
ware (see Section 3.1 in Oil). WFC3/UVIS-G280 spec¬ 
tra span roughly A = 2000 — 6000 A and they have rela¬ 
tively high signal-to-noise ratio (S/N~20) per pixel down to 
A ~ 2000^] (FWHAI ~ 60A at A = 2500A). Uncertainties 
in the wavelength calibration are of the order of 2 pixels, in 
the form of a rigid shift in the pixel space (see Table 1 in 
Oil). _ _ 

As already pointed out by O’Meara et al. (20131 (013 
hereafter), since this sample is based on SDSS optical colour 
selection (JRichards et ah 2002]), it might be biased towards 
bluer colours than a complete sample. However, |Worseck fe| 


Prochaska (2011) have demonstrated that at the redshifts 


of this sample this bias is relatively small (see their Figure 
16). 


Note however that quasar samples constructed from 
UV spectroscopy archives are based both on optical colour- 
selection and additional UV brightness criteria. Typically 
the brightest objects observed in previous HST cycles were 
re-observed when a new UV instrument came online. Ob¬ 
jects observed by the HST Cosmic Origins Spectrograph 
were observed in the post GALEX era, and are thus explic¬ 
itly biased to have bright near-UV and far-LTV magnitudes. 
The selection function of quasars in the UV archives is thus 
unknown and extremely difficult to quantify, but the expec¬ 
tation is that such samples are biased very blue. In contrast, 
the selection criteria for our quasar sample is much cleaner, 
as it is simply an optical apparent magnitude limited sample 
g < 18.5 of quasars at z em ~ 2.4. 


2.1 X-rays 


The X-ray data have been gathered from the ROSAT, 
XMM-iVewjfon, and Chandra archives. We found 2 objects 
(J1253+0516 and J1335+4542) detected in the ROSAT All- 
Sky Survey Faint Source Catalog (RASS-FSC), 6 quasars 
have at least one ROSAT/PSPC ( |White et al. 1994a|b I 
observation, while 2 sources have ROSAT/HRI images 
(J1454+0325 and J1119+1302 with also a PSPC observa¬ 
tions). 

The X-ray fluxes at 0.5-2 keV were calculated from the 
count rates in the ROSAT band (0.1-2.4 keV) and in the 
PSPC band (0.24-2 keV) by utilising a power law spectrum 
with a photon index T = 2.0 modified by Galactic absorp¬ 
tion (|Kalberla et al.||2005|) only. The X-ray properties of 


1 The data generally have S/N exceeding 10 pixel 1 at all wave¬ 
lengths A > 2000A. 


the objects in our sample detected by ROSAT are listed in 
Table IA1I 

Sources without a ROSAT detection might be either 
intrinsically faint, X-ray obscured, and/or highly variable. 
We expect the undetected AGN to have X-ray luminosities 
of the order of 7 x 10 45 3 erg s -1 or lower in the rest-frame 
0.5-2 keV banc0 

Among these 10 quasars, two have additional spec¬ 
tral information from XMM-Newton (J0755+2204 and 
J1119+1302) and one from Chandra (J1220+4608). The soft 
(S = 0.5 — 2 keV) and hard (H = 2 — 10 keV) X-ray lumi¬ 
nosities for J0755+2204 and J1119+1302 are tabulated in 
the SDSS (DR5)/XMM -Newton quasar survey catalog, and 
the values are Ls = 4.60 x 10 44 and Lh = 1.14 x 10 45 erg 
s _1 for J0755+2204, while J1119+1302 has L s = 1.96 x 10 45 
and Lh = 5.00 x 10 45 erg s _1 . 

J1220+4608 was targeted by Chandra with ACIS-S in 
March 2008 with a nominal exposure of 50 ks. The spectrum 
has very low number counts, with a ^ 3<r detection in the 
soft band. In this case, assuming a photon index of 1.8 and 
Galactic Nh, we found that the extrapolated flux in the 2-10 
keV band is 1.36 x 10 -14 erg s _1 cm -2 . This flux corresponds 
to a luminosity of 6.4 x 10 44 erg s~ 4 . 

Summarising, the fraction of detected sources is 21% 
(i.e. 11/53, 10 objects detected by ROSAT/XMM-Newton 
and one additional source with a Chandra observation) with 
a mean soft X-ray luminosity of about 5.6 x 10 45 erg s _1 . 


2.2 Radio 


To estimate the fraction of radio emitters in our sample, we 
matched it with the DR7 quasar property catalog by |Shen| 
et al. (2011). Radio properties are collected from FIRST 


(Becker et al. 19951 and from NRAO/VLA Sky Survey (Con 
don et al.|1998| with a matching radius of 30” (our matches 


are all within 1.3”). We have found 11 quasars in the FIRST 
survey, while one object (J2338+1504) has a detection from 
NRAO, for a total of 12 detections. The radio loudness is es¬ 
timated following Jiang et al. | {2007 ),R = f 6 cm //2500 where 
fscm and /2500 are the flux density at rest-frame 6 cm and 
2500 A, respectively. According to Jiang et al., 11 quasars 
are core-dominant and one is lobe-dominant (.11119+1302, 
radio morphology classification following Jiang et al.J2007). 
Five objects are not in the FIRST footprint (flag = -1 in 
the DR7 catalog, J1259+6720, J1325+6634, J1400+6430, 
J2111+0024, J2136+1029) and they do not have any detec¬ 
tion from NRAO. A summary of the radio properties is given 
in Table A2 Following Kellermann et al. (19891, quasars 


with R > 10 are defined as “radio-loud”. In the WFC3 sam¬ 
ple the radio-loud fraction (RLF) is of the order of 19%. 
This is consistent with the RLF of quasar at similar optical 
magnitudes {] 

In the following, we will present the WFC3 average 


2 The ROSAT RASS flux limit is 5 x 10 -13 erg s _1 cm -2 for 
a mean effective exposure time of 400 sec and T = 2.0, but the 
range of ROSAT exposure times is large, and the sensitivity limit 
is different from field to field ( |Voges et al . 1999[ ). 

3 | Jiang et al. | ( |20(5 t| ) have found that the RLF is about 15% for an 

absolute magnitude at the rest-frame 2500A (M 2500 ) of ~ —28, 
which is consistent with the average M 2500 of the WFC3 sample 
(A /2500 — —29). 
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Figure 1. Mean observed QSO spectrum as estimated from 
a stack of the WFC3 spectra, each normalized to unit flux at 
1450A(see § for details). The black line represents the WFC3 
QSO average spectrum with uncertainties from bootstrap (black 
shaded area). The red line represents the SDSS QSO average 
spectrum (not smoothed to the WFC3 resolution) of our WFC3 
sample with uncertainties. Grey lines represent the single QSO 
spectra, while thin red lines represent SDSS spectra smoothed to 
WFC3 resolution. 


spectrum with and without the 10 objects with R > 10 (a 
stand alone radio-loud quasar stack is not constructed given 
the poor statistics). 


2.3 GALEX 

To extend our wavelength coverage to shorter wavelengths, 
we first considered the GALEX photometry in the DR7 
quasar property catalog. We found that the detection rate 
for the near-UV (NUV at A ef r = 2316A) and far-UV (FUV 
at A e ff = 1539A) bands was 64% (34/53) and 14% (14/53), 
respectively. Given the low detection rate, we have cross 
matched our sample to the GALEX forced photometry cat¬ 
alog (David Schiminovich, private communication), which 
allows us to obtain a detection for almost all sources in 
the WFC3 sample. The GALEX forced photometry was 
not available for only 4 objects (J0751+4245, J1354+5421, 
J1540+4138, J2111+0024), for the rest of the sample we 
have both NUV and FUV bands. The observed GALEX 
forced photometry for the WFC3 sample is listed in Ta¬ 
ble |A3| The average GALEX bands are then corrected for 
IGM absorption and we will outline how these data have 
been corrected for IGM absorption in Section [5.4| 


3 AVERAGE SPECTRUM CONSTRUCTION 

We follow a similar procedure as 013 for the construction 
of the WFC3 average spectrum. Specifically: 

(i) We correct the quasar flux densitjJ^] (/a) for Galactic 


4 In the following we will use the word “flux” to mean the flux 
density (i.e. flux per unit wavelength). 


reddening by adopting the E{B~V ) estimates from Schlegel 
et al. ( 1998| SFD, median reddening value is E(B — V) = 


0.02 mag) and the Galactic extinction curve from Fitzpatrick 


(1999) with Rv = 3.0. The same reddening law has been 
considered to correct the GALEX fluxes. 


(ii) We generate a rest-frame wavelength array with 
fixed dispersion AA. The dispersion value was set to be 
large enough to include at least one entire pixel from 
the WFC3/UVIS-G280 spectra at rest wavelengths A < 
1215A (i.e. AA ~ 6.2A). 

(iii) Each quasar spectrum was shifted to the rest-frame 
and linearly interpolated over the rest-frame wavelength ar¬ 
ray with fixed dispersion 

(iv) We normalized single spectra by their flux at rest 
A = 1450A. 


(v) All the flux values were then averaged to produce the 
stacked spectrum normalized to unity at A = 1450A. 


Recently, Peek & Schiminovich (20131 have found that UV 


colours at high latitudes are best fitted with a |Fitzpatrick 


(|1999[) exti nction curve with Rv — 2.2, while Schlafly & 
Finkbeiner| ( 2011|) found that SFD overestimates redden¬ 


ing by 14% (reddening values should be recalibrated as 
E(B — V) = 0.86 x E(B — V)sfd )• The average WFC3 
spectrum re-estimated following the results outlined above 
is fully consistent with the one constructed with Rv = 3.0 
and E(B - V) from SFD. 


Uncertainties on the observed stack are estimated 
through the bootstrap resampling technique. We created 
5000 random samplings of the 53 spectra with replacement, 
and we applied the same procedure as described above. The 
resulting stack is shown as the solid black line in Figure [T] 
for the full sample, while the resulting uncertainties on the 
stacked spectrum are plotted with a shaded area. The same 
technique is applied to the SDSS quasar spectra and the 
resulting stacked spectrum is presented in Figure [I] with 
the red solid line. For plotting purpose we show the SDSS 
stacked spectrum down to A = 1450A given that the WFC3 
stack already covers wavelengths below this region. 


To bring the WFC3 and the SDSS stacked spectra over 
a common luminosity scale in the vL v plane, we have multi¬ 
plied the final stacks for the mean flux (estimated from the 
WFC3 and SDSS) of the sample at 1450A. For the sake of 
matching the flux scale we have convolved the SDSS spectra 
to WFC3 resolution. The mismatch between the WFC3 and 
SDSS spectra at 1450A is 11% and it is due to variability 
and flux calibration errors in the WFC3 data (see Section 
3.2 in Oil for details). We have thus re-scaled the WFC3 
stacked spectrum to match the (convolved) SDSS at 1450A. 

We have also constructed the quasar average spectrum 
by excluding the 11 objects with R > 10 for completeness. 
The resulting spectrum is shown in Figure [2] where we have 
over plotted the WFC3 average spectrum as a comparison. 
The radio quiet quasar stack is fully consistent with the 
WFC3 stack within the uncertainties. 
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Figure 2. Upper panel: Mean observed QSO spectrum for the 
WFC3 sample (black line) and for the subsample with R < 10 
(red line), each normalized to unit flux at 1450A. Lower Panel: 
Ratio of the R < 10 to the full mean observed QSO spectrum. 



Figure 3. IGM transmission curves as a function of rest-frame 
wavelength. Grey curves represent 10,000 different realisations of 
T\, while the solid and dashed red curves are the stack and la 
dispersion, respectively. The black solid curve is the smoothed 
stack to the WFC3 resolution with the corresponding 1 a disper¬ 
sion (black dashed lines). 


4 IGM TRANSMISSION CORRECTION 


Blueward of Lyman alpha emission in the quasar rest frame, 
absorption from intergalactic HI attenuates the quasar flux, 
both in the Lyman series (creating the so-called HI forest), 


and in the Lyman continuum at rest A < 912 A (e.g. Moller 
fe Jakobsen||1990|). The large abundance of neutral gas at 


z e m < 2.4 is very apparent in our average quasar spectrum 
shown in Fig. [l] With only mild assumptions on the average 
quasar SED, the exponential flux decline at A < 912 A yields 


the mean free path of Lyman limit photons in the IGM 
(013). Here we reverse the question and constrain the quasar 
SED for a range of IGM transmission curves T\ for A < 
1215.67 A. 

For a source at emission redshift z em the effective opti¬ 
cal depth to Hi Lyman series and Lyman continuum pho¬ 
tons at redshift z < z em is determined from the Hi ab¬ 
sorber distribution function in redshift and column den¬ 
sity /(Am, z) = d 2 n/ (dNmdz). The resulting average IGM 
transmission T\ critically depends on the parametrization of 
/(Ah i, z) ( Madau|1995| Meiksin|2006 Inoue et al.|2014 1 and 
is statistical in nature due to the stochasticity of Lyman 


limit systems (Bcrshady et al. 1999 Inoue & Iwata 2008 


Worseck & Prochaska 2011b Quasar composites based on 


low-resolution spectra need to be corrected for Lyman series 
and Lyman continuum absorption of low-column density ab¬ 
sorbers that cannot be identified and corrected by eye. The 
statistical IGM correction strongly depends on the assumed 
absorber distribution parameters and their dependence on 
redshift, which may have resulted in large systematic errors 
in existing quasar composite spectra if the incorrect param¬ 
eters were used or if redshift evolution was not properly 
taken into account. Moreover, the ability to identify weak 
partial Lyman limit systems (T 912 < 0.3) depends on red¬ 
shift, the employed spectra (S/N, spectral resolution, wave¬ 
length coverage, flux calibration), and the intrinsic quasar 
continuum, such that heterogeneous samples are prone to 
ambiguities in the continuum definition and the treatment of 
Lyman limit systems. In particular, we stress that applying 
an average correction for partial Lyman limit systems to de- 
attenuate individual z em < 3 sightlines is incorrect because 


of the stochasticity of Lyman limit absorption! Worseck fc 
Prochaska|2011| ). 


Our large sample at similar z e m ~ 2.4 ensures a suffi¬ 
cient sampling of partial Lyman limit systems, justifying a 
redshift-specific IGM correction function T\ to be applied 
to our stacked spectrum (Fig. [lj. However, the low spectral 
resolution prevents an unambiguous identification of weak 
partial Lyman limit systems in individual spectra without 
knowledge of the underlying quasar continuum. We therefore 
have to correct for the average Lyman series and continuum 
absorption of the whole H 1 absorber population statistically. 
While this is the simplest and most liberal approach, the sig¬ 
nificant flux drop in the Lyman continuum implies a large 
correction to our stacked spectrum, with additional uncer¬ 
tainties related to the parametrization of the IGM. 

In this analysis we consider the recent constraints on the 
z ~ 2.4 IGM presented by Prochaska et al. (2014 and refer¬ 


ences therein) within their range of uncertainty. Prochaska 
|et al. (|2ClG|) introduced a cubic Hermite spline model to 
describe /(Am, 2 )- They performed a Monte Carlo Markov 
Chain (MCMC) analysis of existing constraints on /(Am, z) 
to derive the posterior probability distribution functions of 
seven spline points spaced at irregular logarithmic intervals 
in the range Ahi = 10 12 10 22 cm” 2 . Using the output from 
their MCMC chainsj^J we generated 10,000 realizations of 
/(Ah 1 , z) at z = 2.4 and calculated T\ in the observed wave¬ 
length range with a semi-analytic technique. This modeling 
assumes that the Hi forest is composed of discrete “lines” 


5 Wavelengths are divided by (f + z) to shift the spectra into the 

source rest frame, while fluxes in f\ are multiplied by (1 + z). 6 http://www.arcetri.astro.it/~lusso/Site/Research.html 
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with Doppler parameter b = 24 km s _1 and that the normal¬ 
ization of f(Nn u z) evolves as (1 + z) 0 ' 5 which implies in¬ 
creasing transmission below A = 912 A for the lower-redshift 
Lyman series. Opacity due to metal line transitions was ig¬ 
nored since they contribute negligibly to the total absorp¬ 
tion in the Lyman continuum. In Fig. [3] we plot our suite 
of 10,000 IGM transmission functions with grey lines. We 
then constructed the mean of all of these different realiza¬ 
tions and smoothed this to the WFC3 grism resolution (5 
pixels). Finally, we shifted the observed wavelengths to rest 
frame at z = 2.4 and resampled the transmission functions 
onto the rest-frame wavelength grid of our stacked quasar 
spectrum. This defines our mean IGM transmission func¬ 
tion which is shown as the black curve in Figure [3] Given 
the narrow redshift range of our quasar sample, we did not 
account for redshift evolution in T\. To summarize, by treat¬ 
ing f(N St ,z) within reasonable uncertainties our approach 
yields a range of plausible IGM correction functions, a clear 
advance over previous analyses that assumed a fixed correc¬ 
tion. However, since we account for the total Hi absorber 
population, our IGM correction functions in Fig. [3] assume 
that the column density distribution is well sampled at high 
column densities. 

The procedure to correct the observed WFC3 spectra 
for IGM absorption is outlined as follows: 

(i) We generate a set of 10,000 mock quasar stacks, fol¬ 
lowing the same procedure as in Section [3] by drawing ran¬ 
domly from the 53 quasar spectra to assess sample variance 
allowing for duplications. 

(ii) We then randomly draw one IGM transmission func¬ 
tion from our suite of 10,000. We smoothed this to the WFC3 
grism resolution (5 pixels), and we resampled the transmis¬ 
sion function onto the rest-frame wavelength grid of our 
stacked quasar spectrum. This is repeated for each mock 
quasar stack. 

(iii) We divide the observed spectral flux (/a, obs) by the 
IGM transmission curve 

/A,corr = f\,obs/T\. ( 1 ) 

(iv) The 10,000 mock stacks corrected from IGM absorp¬ 
tion are then averaged to produce the stacked spectrum 
(normalized to unity at A = 1450A). 

(v) The uncertainties on the corrected WFC3 stacked 
spectrum are estimated from the dispersion of these 10,000 
mock stacks. 

The resulting stacks for the full WFC3 and the radio quiet 
samples are shown in Figure [4] with the blue and red lines, 
respectively, and are tabulated in Table[l] The stacked spec¬ 
tra show a softening at wavelengths A < 912A and several 
emission lines. We will discuss these results further in Sec¬ 
tion^ The 1 —a uncertainties on the constructed stack are 
displayed as a light blue shaded area. 

One possible concern about the measure of the uncer¬ 
tainties on the stack is that the bootstrap technique outlined 
above may not have converged far in the blue. In the limit 
where only a handful of the same transmitting spectra actu¬ 
ally contribute to the stack at A < 912A the bootstrap may 
underestimate the true noise in our measurement. In general, 
the overall uncertainty on the stack depends on a variety of 
effects such as the intrinsic fluctuations of the underlying 
quasar continuum (sample variance), our imperfect knowl- 



A rest-frame [A] 


Figure 4. Mean IGM corrected QSO spectrum with uncertainties 
from bootstrap (shaded area) for the full WFC3 (blue) and for the 
radio quiet (red) sample. The black line represents the observed 
WFC3 QSO stack with uncertainties from bootstrap (grey shaded 
area). 


edge of the IGM transmission function (due to uncertainties 
in the /(TVhi, z ); see Figure|4|, spectral noise (i.e. read noise 
and photon counting noise), and most importantly, Poisson 
shot-noise (sample variance) due to the number of LLS and 
partial LLS intercepted by the quasar sightline. 

As we have discussed above, our suite of 10,000 T\ are 
created from an integral over the column density distribu¬ 
tion, which includes uncertainties in the parameters govern¬ 
ing the f(Nm,z), but it does not include the variance due 
to the stochasticity of LLS absorption. 

To this end we conducted a Monte-Carlo simulation 
with a different set of empirical IGM transmission curves 
following the IGM parametrization described in 013. We 
assumed the f(Nni,z) defined in Table 10 in 013 with a 
(1 + z) 1 ' 5 redshift evolution in the normalization (see their 
Equation 6). For each quasar in our sample we drew 100 
sightlines from a parent sample of 5000 sightlines populated 
with 0 < 2 < 3 absorbers and computed mock WFC3 IGM 
transmission spectra in the covered wavelength range, i.e. at 
2000 A< A < 1215.67 (1 + z em ) A. This allows for the small 
redshift evolution in the IGM transmission in the narrow 
redshift range covered by our sample. The result is a set 
of 100 mock spectra per quasar that accounts for the Pois¬ 
son statistics of Lyman limit systems. This method does not 
consider uncertainties in the shape of f(Nm, z), but samples 
/(JVhi,z) properly at high column densities (albeit without 
clustering). The Monte-Carlo simulation is then carried out 
as below: 

(i) We start by assuming the true underlying mean QSO 
spectrum to be the one corrected by the T\ functions. 

(ii) We then create a large number (i.e. 8000) of mock 
samples of 53 quasars by drawing randomly from the 100 
IGM transmission realisations for each quasar in the mock 
sample. 

(iii) We finally calculate the variance of this stack from 
these many ensembles of 53. 
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A rest-frame [A] 


Figure 5. Comparison between the Monte Carlo stack and 1 — 
a uncertainties described in § [4] (red dashed line ) with the one 
estimated by applying the bootstrap technique on the data. The 
red shaded area matches the uncertainties due data bootstrap, 
highlighting the fact that the bootstrapping is converged far in 
the blue. 


This procedure fully encapsulates the stochastic nature of 
IGM absorption, and in particular shot noise due to the 
presence or absence of LLSs. Figure [5] presents the com¬ 
parison between the uncertainties from bootstrap and those 
estimated following the above approach. The amplitude of 
the uncertainty on the Monte Carlo stack matches the one 
estimated by applying the bootstrap technique on the data. 
Our Monte Carlo simulations suggest a ~16% fluctuation in 
samples of 53 due solely to shot-noise in the IGM absorption. 
Since this fluctuation is comparable to the bootstrap error 
in our observations, this strongly suggests that the IGM 
fluctuations dominate the error budget in the stack, and 
that the other sources of variance (i.e. intrinsic variations in 
the quasar continuum and spectra noise) are sub-dominant. 
Most importantly, this comparison also shows that the boot¬ 
strapping is converged far in the blue. In fact, given the nar¬ 
row AGN redshift range, which goes from z m in = 2.282 to 
2 min = 2.599 (mean redshift (z em ) ~ 2.44), almost all spec¬ 
tra contribute appreciably to the total flux at 600 A. We 
note that the stack estimated with this second set of IGM 
transmission curves is in agreement, within the uncertain¬ 
ties, with the one constructed from the 10,000 T\ created 
from an integral over the column density distribution. We 
have thus decided to consider the stack constructed from the 
10,000 IGM transmission curves for the rest of our analysis. 


5 RESULTS 

5.1 Spectral fit and emission lines 

The AGN ionising continuum and its shape in the optical- 
UV is crucial for several reasons, such as interpreting broad 
and narrow emission-lines observed in AGN spectra, and 
their relative ratios. It also defines the BBB (i.e. the bulk of 
the QSO emission), which can be important in interpreting 
the observed relation between optical and soft X-ray fluxes. 


Table 1. WFC3/UVIS Stacked Spectrum corrected for IGM ab¬ 
sorption 


A a 

All 

ft t 

ff (/A,f)° 

RQ 

f\, f 

cr(/A,f ) 

583.351 

2.382 

0.577 

2.405 

0.602 

589.540 

2.247 

0.548 

2.302 

0.580 

595.729 

2.326 

0.542 

2.345 

0.576 

601.919 

2.297 

0.533 

2.248 

0.551 

608.108 

2.234 

0.518 

2.205 

0.548 

614.297 

2.092 

0.492 

2.107 

0.522 

620.486 

2.081 

0.472 

2.049 

0.492 

626.675 

2.090 

0.474 

2.044 

0.495 

632.865 

2.037 

0.462 

1.996 

0.477 

639.054 

1.945 

0.444 

1.911 

0.467 

645.243 

1.937 

0.428 

1.916 

0.448 

651.432 

1.908 

0.418 

1.888 

0.442 

657.621 

1.872 

0.402 

1.864 

0.423 

663.811 

1.866 

0.388 

1.871 

0.413 

670.000 

1.867 

0.388 

1.893 

0.413 

676.189 

1.850 

0.379 

1.890 

0.405 

682.378 

1.848 

0.373 

1.857 

0.394 

688.567 

1.814 

0.363 

1.831 

0.384 

694.757 

1.814 

0.354 

1.847 

0.375 

700.946 

1.811 

0.340 

1.866 

0.370 


Notes. 

a Rest-frame wavelength in Angstrom. 

b Mean IGM corrected flux per A normalized to the flux at 
1450A. 

c Flux uncertainties from our bootstrap analysis (see text). 

(This table is available in its entirety in a machine-readable 
form in the online journal. A portion is shown here for guidance 
regarding its form and content.) 

Table 2. Emission lines properties 


Line 

^obs 

(A) 

Flux a 

(erg s —1 cm -2 A -1 ) 

EW 

(A) 

H I-Ly a+S IV 

938.0 

8.2 ± 2.7 x 10- 16 

0.5 

Lyq+Cm] 

987.6 

4.4 ± 1.1 x 10" 15 

3.0 

Ly/3+Oiv 

1029.7 

9.8 ± 0.4 x 10 -15 

8.5 

Fe Il+Fe III 

1122.0 

1.8 ± 1.0 x 10" 15 

1.3 

Lya 

1216.1 

9.0 ± 0.9 x 10" 14 

74.0 

Si iv+O iv] 

1397.5 

6.7 ± 1.4 x 10- 15 

7.2 

Civ 

1544.6 

1.7 ±0.2 x 10" 14 

17.8 

Hen 

1635.7 

1.2 ± 0.3 x 10" 15 

2.0 

O hi] 

1663.7 

7.3 ± 1.0 x lO" 16 

1.2 

Aim 

1861.0 

1.8 ±0.5 x 10" 15 

2.6 

Cm] 

1907.3 

9.7 ± 0.9 x 10- 15 

15.5 


Notes. 

a Normalized fluxes are multiply by the average flux of the 
WFC3 sample at 1450A. 
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Figure 6. Left panel: Average WFC3 and SDSS spectra normalized at 1450A with fits to power-law continuum (dot-dashed line). Five 
continuum windows (from T02) are shown as black horizontal lines below the spectrum. The power-law continuum fit exhibits a break 
at A ~ 920A, with a flatter (softer) spectrum at shorter wavelength. Right panel: Zoom in of the EUV region (600—912A). The dashed 
line is the power-law continuum by employing the average EUV spectral slope. 


In order to infer the continuum slope of the WFC3 
stack, we have to avoid contamination of broad emission 
lines, broad absorption features and extended wings of emis¬ 
sion line^] especially at A > 1216A where the emission lines 
are more prominent. 

We have thus fitted the UV stack in each wavelength 
window free of strong features by means of a single power 
law. We adopted the same intervals as T02, but given that 
the uncertainties dramatically increase at A < 912A we de¬ 
cided to restrict the fit of the continuum redwards of 912A 
(i.e. 1095-1110, 1135-1150, 1450-1470, 1975-2010, and 
2150—2200 A). We note that the continuum window blue- 
ward of Lyct (i.e. 1095— lllOA) may be contaminated by 
the Fell+Fem emission line. Given our limited resolution 
we decided to keep it for consistency with T02. Results are 
unaffected if we neglect this window. We have considered as 
the error for the \ 2 minimization procedure the uncertain¬ 
ties in the stack, which are rather uniform at A > 912 A. 
Figure [6] shows the rest-frame stacked spectrum extending 
from 600 A to 2500 A and the power-law fit to the contin¬ 
uum of the form f\ oc A“*, where the best-fit power law 
index is a v = —0.61 ± 0.0f[^] (dashed line). 

Our best-fit slope is consistent with the one estimated 
by T02 in the FUV (—0.69 ± 0.06 at A > 1216A, see their 
Table 1), and by S12 (—0.68 ± 0.14) in the same region (see 
Figure |8|. We note that the continuum regions adopted by 
S12 are the same as T02. Vanden Berk et al. ( 2001| ) also 
measured the continuum slope of a quasar composite utilis¬ 
ing 2200 spectra from SDSS, which span a redshift range of 
0.044 ^ 2 ^ 4.789. They find a slope of a v = —0.44 ± 0.10 


7 Emission from host-galaxy star light is not considered in the 
fit since our stack covers the rest-frame UV range of very bright 
AGN, where the expected contribution of the galaxy should be 
negligible. 

8 In the following we will refer to a„ only. The relation between 
the fluxes in wavelength, f\ oc A'**, and frequencies, f u oc r” 1 ', 
is a„ = —(2 + a x ). 


(no IGM absorption correction is applied), which is signifi¬ 
cantly different than the one we measured at 2.4<r. The dif¬ 
ference can be attributed to the different continuum regions 
adopted by Vanden Berk et al. (20011 (i.e. 1350—1365 and 
4200— 4230Ay 

We have also applied different continuum regions to 
gauge the dependence of the fit on the windows adopted. We 
considered the intervals listed by Decarli et al. (2010) for a 
sample of 96 quasars at redshift lower than 3 with spectra 
collected from several instrument (e.g., SDSS, ESO/NTT, 
Nordic Optical Telescope, and HST). Specifically, the al¬ 
ternative windows adopted are 1351 — 1362, 1452—1480, 
1680-1710, 1796-1834, 1970-2010, and 2188-2243A. To 


also fit the UV part of the spectrum, we add to these the 5 in¬ 
tervals at A < 1200A by T02. We find a slope of — 0.63±0.01. 

From Fig[6] it is apparent that a single power-law pro¬ 
vides an excellent fit up to ~ 900A where the continuum 
exhibits a break. A flatter (softer) spectrum is present at 
shorter wavelength, which is not modelled by the simple 
single power-law. Solely for the purposes of a quantitative 
estimate of the position of the break, we have modelled the 
stacked spectrum by adding an exponential attenuation to 
the power-law as 


, JA“\ if A > Ai, 

/A0C |A“a x e- (Ai >- A ) /A G if 600A < A < Aj, 

where At, and A/ represent the break and the attenuation 
wavelengths, respectively. The best-fit power law index is 
a v = —0.65 ± 0.01, the break occurs at Aj = 922 ± 31 with 
an attenuation factor A/ ~ 450 not well constrained given 
the large uncertainties on the stacked spectrum in the far 
blue. 

The break measured from our stacked spectrum is at 
bluer wavelengths with respect to the ones estimated from 
composites at much lower redshift. S12 found the break 
at A ~ lOOOA, whilst T02 located the wavelength break 
around 1200—1300A(a comparison among the various com¬ 
posites is given in Figure [8]). The difference from the latter 
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Figure 7. Histogram of the EUV slopes estimated from the pro¬ 
cedure outlined in § |5.1| The solid and the dashed lines represent 
the mean and 1 — a dispersion, respectively, while the dot-dashed 
line denotes the median. 


is due to the poor IGM correction at those wavelengths, 
which artificially produces a break around the Ly a. We will 
analyse this issue further in Section |5.2| Additionally, the 
composite shown by S04 is much steeper than what we ob¬ 
served in the same wavelength range (i.e. a„ = 0.56to'28 
A = 650 — 1150A). 

Given that emission lines are much fainter in the EUV 
region than in the FUV/NUV and that uncertainties far in 
the blue are significant, we have considered the full stacked 
spectrum (continuum + lines) at A < 912A and found that 
the EUV slope is —1.65 ±0.07. The error on the EUV spec¬ 
tral index from the fit of the average spectrum is quite 
small given the uncertainties on the WFC3 stack itself at 
these wavelengths. Additionally, the flux in the EUV sys¬ 
tematically depends on the IGM correction applied, such 
that fluxes and errors are strongly correlated. The simple 
X 2 fitting technique is not appropriate in the EUV, and the 
bootstrap procedure is required. We have thus estimated the 
EUV slope from each realization in the bootstrap described 
in Section [4] Figure [7] shows the distribution of the EUV 
slopes of mock stacks. We found qeuv = —1.70 ±0.61 (con¬ 
tinuum + lines), which is a more reasonable value of the ac¬ 
tual uncertainty. We will quote this as the final EUV slope. 
Lastly, we comment that it is evident from the right panel 
of Fig.[6]that a single power law does not seem to be a satis¬ 
factory description of the region below 912A. For example, 
there is no good reason why the feature at ~ 730A should be 
considered an absorption feature. However, given the large 
uncertainties of our stack and to compare oeuv with pre¬ 
vious evaluations from the literature, we refrain either to 
employ more complicated functions or to define continuum 
windows to fit this region. We further discuss this issue in 
Sect. 15.61 

We also identify most of the emission lines in our 
stacked spectrum usually seen in optical spectra of high- 
redshift QSOs such as Ly/3, Lya, Siiv+Oiv], Civ, and the 
semi-forbidden line of C ill] 1909. A number of weak lines 
show up in the average spectrum at A < 1216A, including 
NeVm+Olv 772, Om 831, Lyy+Cm] 873, and Fell+Fem 


1123. Blended lines from high-ionisation states such as Oiv 
608, Ov 630, Nhi 685, and Om 702, which are important 
diagnostics for investigating the physical conditions of broad 
emission line regions, may also be present, although it is im¬ 
possible to reliably measure their strengths given the noise 
in our stacked spectrum at blue wavelengths. Weaker lines 
in the FUV include Hen 1640, and Om] 1663. The prop¬ 
erties of the emission lines for which the uncertainties on 
the stacked spectrum are small (i.e. A > 912A) are listed in 
Tabled 

Summarising, we estimated the continuum slope in the 
NUV/FUV {a v = -0.61 ± 0.01) and in the EUV (a„ = 
— 1.70 ± 0.61). We confirm the presence of a break in the 
quasar stacked spectrum by employing a completely inde¬ 
pendent sample at high redshift. We also emphasise that 
previous estimates of the break were performed with very 
few spectra (~ 10 in T02 and 3 — 12 in S12) contributing at 
short wavelengths, whilst all 53 spectra in our high-redshift 
sample are contributing at 600A. 


5.2 Comparison with Previous Quasar 
Composites 

In this section we will perform a detailed comparison be¬ 
tween our z em ~ 2.4 quasar average spectrum with previous 
works in literature. We caution the reader that a comprehen¬ 
sive and consistent explanation for all the differences hinges 
on multiple factors as, for example, sample selection biases. 

In Figure [8] we compared our WFC3 average spectrum 
with the AGN composites found from SDSS by Vanden Berk 
et al.| ( |2001[ ) , from HST by T02, from COS by S12, and 
from FUSE by S04. The Vanden Berk et al. (20011, T02, 
and S12 composites are normalized at 1450A, whilst the 
S04 composite is, instead, normalized to our IGM corrected 
mean stack at A = 1114A (logiz = 15.43). We caution that 
the latter normalization is more problematic, because it is 
subject to uncertainties in the Lya forest IGM correction 
employed by S04. 

Several important points emerge from the comparison 
in Figure [8] First, it is apparent that the IGM correc¬ 
tion employed by T02 is significantly underestimated at 
A < 1216A. In general, T02 attempted to identify individual 
LLSs by eye and used estimates of the Lyman limit optical 
depth to correct spectra with strong absorption. Whereas for 
lower column density systems, a single statistical correction 
for unidentified Lyman limit absorbers in the Lya valley 
( Moller fc Jakobsen| 1990) was applied per spectra using the 
column density distribution 


d 2 n 

dzdN 


oc (1 + z) 1 N 


(3) 


where n is the number of lines, a is the redshift, and N is the 
column density of neutral hydrogen. Values for the /3 factor 
considered by T02 are fi — 1.83 for 3 x 10 14 < IVhi < 10 16 


cm 2 and /3 = 1.32 for Nui > 10 16 cm 2 (Petitjean et al. 
1993|. For 3x 10 12 < Nh, < 3xl0 14 cm" 2 T02 used /3 = 1.46 


(Hu et al. 1995). 

The application of this formula leads to a correction 
characterized by a stepwise behaviour, in which the opacity 
decreases toward shorter wavelengths (as already discussed 


in Binette et al. 20031. We computed the absolute i-band 


© 2002 RAS, MNRAS 000, |T|fT9l 



























10 


E. Lusso et al. 


A [A] 

1000 


Telfer+02 
Shull+12 
Scott+04 
Vanden Berk+01 
Lusso+15 


lSj 



V [Hz] 


Figure 8. Mean IGM corrected QSO spectra with uncertainties (keys as in Fig. [I]and[4| in the rest-frame log v — vL v plane. Rest-frame 
wavelengths (in Angstrom) are plotted on the top x-axis. The T02, S12 and Vanden Berk et al. (2001) composites are normalized to our 
average spectrum at 1450A, and are shown for comparison with the orange, green, and cyan solid lines, respectively. The S04 quasar 
composite (light red solid line) is normalized to our stack at A = 1114A (log v = 15.43). 


magnitude (Mi(z = 2Q for the T02 sample and these values 
are plotted as a function of redshift in Figure [9] There are 
about 20 quasars at 2 > 2 (contributing to the flux at A ~ 
300A in the SED analysis) in the T02 sample and for these 
objects their correction is lower by a factor of ~ 2 in flux at 
A = 600A with respect to ours (see Fig. 3 in Binette et al. 
2003). 

The spectral slope measured by S04 in the EUV is 
a v = —0.56±o;||, which is at variance with the spectral 
shape shown by our WFC3 stacked spectrum (a ~ —1.7 and 
with the other spectra from different samples) in the same 
wavelength range. Like T02, S04 adopted a similar proce¬ 
dure to correct for IGM absorption. LLSs were individually 
identified by eye and a single statistical correction for the 
line-of-sight absorption due to the Lya forest and the Lyman 
valley was applied on single spectra as in Eq. <§■ At variance 
with T02, the parameters to correct for the Lya forest ab¬ 
sorption are j3 = 2 for 1.6 x 10 12 < JVhi < 2.5 x 10 14 cm~ 2 , 
and j3 = 1.35 for 2.5 x 10 14 < A r ni < 5 x 10 16 cm -2 . For the 
redshift distribution parameter S04 used 7 = 0.15. This dif¬ 
ferent parametrization leads to an IGM correction that goes 
from 5% to 10% in the wavelength range 1200—600A at 
2 ~ 0.16 (see their Fig. 4). Any judgment about this IGM 


9 The absolute i-band magnitudes normalized at z = 2, K- 
corrected following Richards ct al. 2006b for the WFC3 sam¬ 
ple have been taken from the DR7 quasar catalog by Shen et al. 
( 2011 ). 


correction being underestimated is not trivial given the lower 
redshifts of the S04 AGN sample, compared with those of the 
WFC3 sample as shown in Figure [9] However, biases in the 
S04 sample, and subtleties in correcting for LLSs may con¬ 
tribute in the differences between the FUSE and the WFC3 
spectras. 

We have also compared our WFC3 average spectrum 
with the one by S12. The COS composite does not signif¬ 
icantly differ from our spectrum within the uncertainties, 
but this similarity cannot be easily explained given the dif¬ 
ferent average redshift, the sample selectioiF] and the IGM 
correction employed by S12, which is solely performed by 
visual spectral inspection. The advantage of the S12 sample 
is the high S/N and the high spectral resolution, which al¬ 
low them to fit the local continua of their objects (correcting 
for identify LLSs). Although S12 missed some partial LLSs 
with small Lyman continuum optical depths, the composite 
constructed by S14 does not significantly differ from the one 
in S12. This is due to the fact that the partial LLSs that 
S14 found have very low columns, logA^Hi < 16 (mostly 
in the 15.0-15.5 range) and their Lyman continuum opacity 
was negligible. Additionally, S14 used the pattern of higher- 
Lyman series lines rather than the Lyman edge. The fitted 
EUV composites of S12 and S14 are in fact quite similar: 
(a„) = -1.41 ± 0.21 in S12, while S14 found -1.41 ± 0.15. 

10 The COS sample was selected from the archive as the best 
available high-S/N spectra of AGNs in January 2011. 
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Figure 9. Absolute i-band magnitude (normalized at z = 2, K- 
corrected following |Richards et al.|2006b[ > as a function of redshift. 
Shaded areas indicate the redshift and magnitude ranges for the 
different samples, estimated from the 16 th and 84 th percentiles. 
Large filled circles represent the median for the different samples: 
Our WFC3 sample (black), T02 (orange), S04 (red), S12 (green), 
and S14 (magenta). Indicative values of the black hole masses (in 
units of Mq) are plotted on the y-axis on the right. Mbh is esti¬ 
mated via Agdd = ^boi/^Edd assuming the average Agdd — 0.35 
for the WFC3 sample. We have estimated the relation between 
Lboi and Mj(z = 2) to be logLboi = —10.03Mi(z = 2)/26 + 36.32 
by fitting the sources in the DR7 quasar catalog. 


Table 3. EW comparison 


Line ionisation EW(line) / EW(Lya) 

energy (eV) a ion 

-1.2 -1.6 -2.0 -1.7 (Lusso+15) 


Lya 

13.6 

1 . 

1 . 

1 . 

1 . 

Cm] 

24.0 

0.1 

0.10 

0.13 

0.21 

Hen 

54.4 

0.09 

0.06 

0.03 

0.03 

Civ 

48.0 

0.78 

0.44 

0.25 

0.24 


Figure 10. Comparison between the three types of SEDs adopted 
by |Baskin et a l. (2014J with our WFC3 average spectrum. The 
green, red, and magenta solid lines correspond to Uion = —1-2, 
-1.6 and -2.0, respectively. The blue points represent the IGM 
corrected mean GALEX forced photometry, where horizontal bars 
indicate the GALEX band-passes. 

1 keV (see their Section 2.2). These three SEDs are plot¬ 
ted in Figure |10| We compare the observed and predicted 
line ratios, rather than the equivalent widths (EW) of single 
lines, since the average covering factor of the BLR in the 
WFC3 sample could be different than the covering factor 
of 0.3 assumed by Baskin et al. In order to calculate the 
predicted line ratios, we use the Baskin et al. model with 
solar metallicity, and a BLR distance where the line emis- 
sivity peaks (see their Figure 5), which is roughly the ex¬ 
pected EW if one assumes a BLR which spans a range of dis¬ 
tances. The predicted Cm]/Lya, Hell/Lya, and Civ/Lya 
for aw = —1.2, —1.6 and —2.0 are listed in Table[3] together 
with the observed line ratios. The observed Hell/Lya and 
Crv/Lya are consistent with the dion = —2.0 model, which 
is softer than our derived slope of-1.7, though within the un¬ 
certainties. The observed Cm]/Ly a is a factor of two larger 
than expected for all assumed values of cuon- 


We also note that, in the S14 analysis only ~ 20 spectra 
contribute at 700A, while few of them do not cover 912A 
such that the continuum normalization may be problem¬ 
atic. Given the bias in the UV samples, and the different 
approach in the IGM absorption correction, it is unclear 
what significance to attach to this agreement. 


5.3 Comparison with photoionisation models 


In this Section we compare the line ratios observed in the 
WFC3 average spectrum and the line ratios predicted by 


photoionisation models. Baskin et al. (2014) presented a 


radiation-pressure-dominated hydrostatic solution for gas 
in the broad line region, using the photoionisation code 
CLOUDY (Ferland et al.|1998 l We defer the reader to that 
paper for details of the physical model. Here, we compare 
their predicted line ratios with our estimates, for several 
prominent BLR emission lines. The predicted line ratios are 
based on three types of SEDs, differing from each other in 
the ionising slope ai on at energies between 1 Rydberg and 


5.4 Comparison with GALEX 


We have also increased our coverage at shorter wavelengths 
by considering the GALEX photometry of the WFC3 sample 
(see Section 2.31. These data have been corrected from IGM 
absorption as follows 


(i) We generate a set of 6000 mock quasar samples, follow¬ 
ing a similar procedure as in Section[3] Each sample contains 
53 NUV and FUV fluxes (we again allow for duplications), 
and we compute the mean of these NUV and FUV fluxes. 

(ii) We normalize the average GALEX fluxes to rest frame 
1450A (as for the WFC3 spectra). The 1450A flux is es¬ 
timated from each average WFC3 spectrum of the mock 
quasar samples. 

(iii) We then randomly draw one IGM transmission func¬ 
tion from our suite of 10,000 and we integrate this function 
over the GALEX filter curves. This is repeated for each mock 
quasar sample. 

(iv) The 6000 mock stacks are thus corrected from IGM 
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Table 4. Accretion disc model parameters 


Ref. model 

M B h 

A E dd 

^corona 

a 


(Mq) 


(R s ) 


ADI (blue) 

6 x 10 9 

0.35 

8 

0.8 

AD2 (red) 

3 x 10 9 

0.70 

20 

0.3 

AD3 (green) 

1.2 x 10 10 

0.17 

8 

1.0 

AD4 (black) 

6 x 10 9 

0.35 

20 

0.0 


Notes. 


All models have the following parameters fixed: T = 1.9, 
kT e = 0.15 keV, f p \ = 0.2, r = 18, and r ou t = 10 5 f? g . 


absorption and averaged to produce the final mean NUV 
and FUV fluxes. 

(v) The GALEX fluxes are then multiply by the average 
flux at 1450A of the WFC3 sample. 

(vi) The uncertainties on the corrected GALEX fluxes are 
estimated from the dispersion of these 6000 mock stacks. 

Both IGM absorption corrected NUV and FUV bands are 
plotted in Figure [To] We ignore quasar variability when 
stacking and plotting the results together with the stacked 
spectra. The GALEX NUV band almost covers the WFC3 
spectra in the Lyman continuum, hence the fact that the 
NUV flux is consistent with the WFC3 stack within the un¬ 
certainties further check our results. 

Nonetheless, the error bars on the NUV and FUV fluxes 
are significant and do not allow us to further constrain any of 
the three SEDs, even though the FUV flux seems to exclude 
the SED with ai on = —1.2 consistently with the comparison 
with the predicted EWs. 


5.5 AGN Accretion disc models 


A basic prediction of simple accretion disc models is that the 
disc temperature decreases as the black hole mass increases 
( Shakura fe Sunyaev||1973 1, thus 


/ GMM\ 1/4 
\ 4nar 3 J 


6.3x 10 5 


/ M 
VAfs 


1/4 

m 8 - 1/4 



(4) 

where M/Me is the Eddington ratio (A E dd, the accretion 
rate normalized to the Eddington accretion rate), Mg = 
Mbh/W 8 Mq, and R g is the gravitational radius ( R g = 
GM/c 2 ). For a scale r of 6 R g and an Eddington ratio 
A Ed d = 0.1, the disc temperature goes from ~ 5 x 10 5 K to 
~ 8.7 x 10 4 K for a Mbh of 10 6 and 10 9 Mq, respectively]^] 
The disc temperature sets the peak of the BBB and, thus, 
we expect to see the location of the break changing as a 
function of Mbh- 

The UV composites by T02 and S12 show a break be¬ 
tween Lya and 1000 A, whilst the quasar composite by S04 
does not show any. This might be consistent with the fact 
that the S04 is the lowest luminosity/redshift sample and 


11 Equation ( 4 ] is valid in the Newtonian limit, and therefore 
may be not accurate to estimate the temperature in the disc, 
which is highly relativistic. See Eqs. (3)-(5) and Table 1 in Laor| 


& Davis 


I 2011 1 for a more accurate approach. 


presumably it has a lower Mbh on average than the S 12 and 
T02 samples. Unfortunately, black hole mass values are not 
available for most sources in the T02, S04 and S12 samples. 
S04 have compiled black hole mass values from the liter¬ 
ature for 22 objects in their sample (with a median Mbh 
of 1.6 x 10 s Mq) and they found a significant correlation 
between the spectral index and black hole mass, with the 
spectral slope being softer for higher Mbh, which may run 
in the expected direction. However, as we have discussed in 
the previous Sections, both the presence and the position 
of the break depend strongly on the IGM correction consid¬ 
ered. 

Furthermore, the standard black body disc model de¬ 
picted above, does not reproduce the observed soft X-ray 
emission seen in AGN, consequently other physical param¬ 
eters are required such as black hole spin and/or substan¬ 
tial Compton upscattering (Laor et al. 1997| but see also 
|Capellupo et al.||2015| ). 

We further analysed our WFC3 average spectrum in 
the context of accretion disc models. In particular we have 
considered the publicly available energetically self-consistent 
model, OPTXAGNF, developed by Done et al. (2012) within 
the XSPEC spectral fitting package. Their model contains 
three distinct spectral components, all powered by the en¬ 
ergy released by a single accretion flow of constant mass ac¬ 
cretion rate, M, onto Mbh: an outer and inner disc, and an 
X-ray corona. The outer disc emits as a (colour temperature 
corrected) blackbody, the inner disc is where a fraction of 
the disc emission is Compton upscattered, while the X-ray 
corona is where also a fraction of the emission is Compton 
upscattered producing the power-law tail at high energies 
(see Section 4 in |Done et al.||2012| for further details). This 
model is set by 9 parameters: Mbh, A E dd, black hole spin 
(a), radius of the X-ray corona (r cor ona), outer accretion 
disc radius (r ou t), electron temperature ( kT e ), optical depth 
(r) of the Comptonization region, photon index (F), and 
the fraction of energy dissipated as a hard X-ray power law 
from v corona to the innermost stable circular orbit ri sco (/ p i). 

Our data deliver constraints on Mbh and A E dd, which 
have been collected from the DR7 quasar catalog (Siren et al. 


2011) for all quasars in the WFC3 sample. The average Mbh 
estimated from the Civ line is 6 x 10 9 Mq, whilst the Ed¬ 
dington ratio is 0.35. The reliability of utilising the C IV 
line is controversial, since C IV can be severely affected by 
non-virial motions such as outflows, winds, and strong ab¬ 
sorption (e.g. |Baskin &; Laor||2005| |Shen et al.||2008| ). We 
will refer to this model as ADI. We have thus investigated 
two additional models with half (AD2) and double (AD3) 
the average A/bh, where the average A E dd has been re-scaled 
consequently. The photon index has been fixed to 1.9, which 
is usually observed in quasars (e.g. Piconcelli et al. 20051. 
For the other parameters we do not have any constraint from 
the data, and hence we have to assume values in order to re¬ 
produce both NUV/FUV/EUV and X-ray in a reasonable 
way. We have thus fixed the electron temperature kT e at 
0.15 keV (Gierlinski & Done 20041, f p i= 0.2, r =18, and 
the default outer accretion disc radius, r out = 10 5 _R g . The 
parameter values for the reference models are given in Ta- 
ble|4] and plotted in Figure [ll| All models are normalised to 
1450A. We note that all the models are only meant to pro¬ 
vide a reasonable description of the observables, since we 
did not attempt any simultaneous fit of the optical-UV/X- 
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Figure 11. Accretion disc models. Keys as in Figure [To| The blue (ADI), red (AD2), and green (AD3) dashed lines represent the three 
reference accretion disc models for Mbh = 6 X 10® with Andd = 0-35 (a = 0.8), 3 X 10®M@ with Andd = 0.70 (a = 0.3), and 1.2 X 10 10 
with AEdd = 0.17 (a = 1.0), respectively. The black dashed line is the same as ADI but with a = 0 and r co rona = 20 (AD4). All models 
are normalised to 1450A. X-ray data for the 10 quasars detected by ROSAT are plotted with filled black points. The XMM-Newton soft 
and hard luminosities for J0755+2204 and for J1119+1302 are plotted as black diamonds. The black arrow represent the ROSAT flux 
limit at 1 keV. 


ray data. X-ray data for the ten quasars detected by ROSAT 
and XMM-Newton are shown in Figure [Til for completeness. 

A key aspect of these models is that the peak of the 
BBB is set by the black hole mass, spin and mass accretion 
rate through the outer disc. Interestingly, the AD4 model 
(i.e. average Mbh and AEdd and spin zero) do not provide a 
good representation of the FUV/EUV region of our average 
spec trunp*] w ith the BBB being colder than observed (see 
also [Lawrence 20121. This problem can be solved if we in¬ 


stead consider the AD2 but with the BH spin set to zero. If 
the BH spin increases, the distance of the innermost stable 
circular orbit reduces, and this effects the disc temperature 
as described by Eq[4] 

Summarising, if we consider the BH masses from C IV 
to be correct, we need to have non zero values of the BH spin 
to match the data, with higher values of the BH spin as the 
Mbh increases (see Trakhtenbrot 2014 for similar results, 
i.e. black hole masses above ~ 3 x 10 9 Mq at 2 ~ 2 — 3 are 
consistent with very high spins). We do not need BH spin if 
we instead consider a factor of 2 mass lower, on average. The 
latter is in agreement with previous results that have found 
BH mass values from Civ systematically higher, once cali¬ 
brated with virial mass estimators based on the Mgn (e.g. 


12 The radius of the corona has been fixed to 20 R g , instead of 
8 as in ADI. 


Shen et al.|2008[|Park et al.|201 3] but see also ITrakhtenbrot] 
fe Netzer|2012 1. 

In principle, the combination of reliable Mbh and 
FUV/EUV spectra (properly corrected for IGM absorption) 
may provide hints on the quasar BH spin. 

It is not possible to test if the location of the break 
changes as a function of Mbh within the WFC3 since it 
has a narrow luminosity range. We thus need to have a 
comparison sample at low redshift/luminosity, hence lower 
Mbh- The FUSE composite by S04 is constructed with low 
redshifts/luminosities AGN (see Fig[9| and interestingly it 
shows a hard spectral slope in the EUV with no break con¬ 
sistent with the BBB being produced by small Mbh- This 
might point towards the direction of the predicted trend, 
but further studies are needed. In particular, a natural ex¬ 
tension of the present analysis is to cover bluer wavelengths 
with high resolution quasar spectra. 


5.6 Quasar emissivity and photoionisation rate 


The UV background is governed by the ionising emissiv- 
ities of the relevant source populations (typically quasars 
and star-forming galaxies) and cosmological radiative trans¬ 
fer in the IGM (e.g. |Faucher-Giguere et al.|[2009[ |Haardt fc| 


Madau 20121. The specific comoving emissivity of the quasar 
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Figure 12. As Fig. [ 8 ] but including recent broken power law 
parameterizations of the quasar continuum (Cowie et al. 1120091 

| Faucher-Giguere et al.||2009| |Haardt &; Madau|2012[ >. The SEDs 

have be en norm ali ze d at the adopted rest frame wavelengths 
(Faucher-Giguere et al. 2009 4400 A) or assumed at 1450 A if 


power law 


not given by the authors ( |Cowie et al. 2009 H aardt &; Madau| 

|2012[ ). All parameterizations underestimate the quasar Lyman 
limit flux of our stacked spectrum and recent quasar composite 
spectra (S12, S14). 


population at redshift z and frequency v, 

poo 

tv(v,z)= / (j>(L,z)L u (L,v)dL , (5) 

critically depends on the adopted quasar luminosity function 
4> (L, z), its faint end (L m i n ), and the quasar SED ( L , i/)). 
In particular, since the quasar luminosity function is typi¬ 
cally determined in the NUV-optical, the specific emissivity 
at the Lyman limit e„,gi2 depends on the chosen FUV and 
EUV SED parameterization, which is commonly taken as a 
luminosity-independent broken power law based on quasar 
composite spectra (Z97, T02). Resulting uncertainties in the 
specific emissivity are rarely discussed (see Faucher-Giguere 
et al.|2008 for an exception), or limited to variations in the 


EUV spectral index ( [Faucher-Giguere et al.| |2009). Conse¬ 
quently, current estimates of the specific Lyman limit emis¬ 
sivity of quasars still vary by a factor ~ 2 at z = 3-4 (Hop- 
kins et al.|2007[|Siana et al.|2008||Cowie et al.|2009[|Masters 


et al. 20121. 


Figure [12] compares our stacked spectrum and various 
quasar composites to broken power-law parameterizations 
used to estimate the specific quasar Lyman limit emissiv¬ 
ity and/or the quasar contribution to the UV background 
(Cowie et al. 2009 Faucher-Giguere et al.||2009 Haardt & 


|Madau||2012| ). Most previous SED parameterizations have 
underestimated the quasar Lyman continuum flux because 
they have been based on the T02 composite spectrum with- 



as correction factors to the quasar Lyman limit emissivity 


£1^,912(2) for a fixed luminosity functiorp^] We point out that 
Cowie et al. (2009) computed a power-law slope of -1.35 at 
700 < A < 2300 A obtained from GALEX photometry of a 
z e m ~ 1 AGN sample. This implies a 700A to 2300A flux ra¬ 
tio of 0.20. Cowie et al. then assumed the actual SED below 
the break to be /„ oc ix~ 4 , which returns an extrapolated 


flux ratio /A912/ fv, 2300 = 0.26. In Fig. 12 we have taken the 
break at exactly 912^p*| We also note here that the mean 
2300A luminosity for the objects in the Cowie et al. sample 
is much lower ( vL u ~ 5 x 10 44 erg s _1 ) than the one of our 
WFC3 sample. Adopting a standard concave SED shape 
instead, their e„,gi2 values increase by a factor of ~ 2. While 
this correction makes their values consistent with higher es¬ 


timates of £^,915 

(Hopkins et al. 2007 Siana et al. 

2008 

Mas- 

ters et al.|2012 

, we note that Siana et al. (2008 

1 and 

Mas- 


ters et al. ([2012 I assume a low flux ratio /A912/./A1450 = 0.58 
based on the T02 composite to convert their £^1450 to £^,912. 
Our SED yields a 30% increase to their estimates, such that 


£^,912 ~ 5 x 10 24 ergs 1 


Hz 1 Mpc 3 at z ~ 3.2. With this 


correction, quasars may account for the total specific Lyman 
limit emissivity estimated from the Lya forest at z < 3.2 


(Becker & Bolton 2013), although a comparable contribu¬ 


tion from star-forming galaxies is allowed within the current 
large uncertainties (Nestor et al.|20l3 Mostardi et al.|2013 


Becker & Bolton 20131. The quasar contribution seems to 
decrease substantially to z ~ 4, as the low (corrected) value 


of e„,9i2 ~ 1.8 x 10 24 ergs 1 Hz 1 Mpc 3 from Masters et al. 


(2012) is just ~ 20% of the total value inferred from the Lya 
forest (Becker & Bolton 2013, although see Glikman et al. 
2011 for a ~ 15 times higher estimate of £^,1450 at z ~ 4). 


From the Lyman limit emissivity we can estimate the 
quasar contribution to the UV background photoionisation 
rate 


Fhi(z) = f 


47 T (u, z) 


hv 


<7hi [y) d v 


( 6 ) 


where J v (y, z) is the mean intensity of the ionising back¬ 
ground at redshift 2 and ohi (u) ~ 0912 (ix/ix 9 i2)~ 3 is the Hi 
photoionisation cross-section with (7912 = 6.33 x 10~ 18 cm 2 
at the Lyman limit frequency 1x912. At high redshifts 2 > 4 
the mean free path of Lyman continuum photons in the 
IGM, Amfp (ix, 2), is much shorter than the Horizon, such 
that the UV radiation field at any given point is dominated 
by local sources within a sphere of radius A m f p , and cosmo¬ 
logical redshift effects in J„ (ix, 2) can be neglected. In this 
‘local-source approximation’ the intensity of the UV back¬ 
ground is 


1 3 

J v (ix, 2) ~ —A m fp (ix, 2) (1 + 2) e v {y,z) (7) 

47r 


(Madau et al. 1999| Schirber fe Bullockj 2003 Faucher 


|Giguere et al.||2008| |Becker fe Bolton||2013| , with A m f p (e„) 
given in proper (comoving) units. The proper mean free path 
to Lyman limit photons is well described by a power law 


Amfp ,912 = 37 


1 + 2 


Mpc 


( 8 ) 


13 Note that |Haardt Sc Madau| ( |2012[ > simplistically assume 
e z/, 912 ( 2 ) and L u (iy) to be independent. 

14 The broken power-law quasar SED is then f u oc v ~ 1,45 (f u oc 
v- 1 ) at A > 912 A (A < 912 A) 
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at 2.3 < 2 < 5.5 (Worseck et al. 20141, whereas its frequency 
dependence is 

/ \ 1.5 


Amfp (^j z) — A m fp,912 


\^912 / 


(9) 


for a column density distribution /(Am, 2 ) oc N B 


I Faucher-Giguere et al. 


2008 


With our power-law param¬ 


eterization of the quasar SEE^_] /„ oc v aEVV at A > 912 A, 


we can write the comoving specific emissivity of quasars as 

Q EUV 

( 10 ) 


( X V V 

( U,Z ) = e„,9i2 ( 2 ) - 

V ^912 / 


Combining Equations 0 ( |10[ ) and integrating Equation 
the quasar contribution to the UV background photoionisa¬ 
tion rate is 


r H , ~ 4.6 x 10 -13 s _1 


Ci/,912 


10 24 ergs -1 Hz 1 Mpc 


1+2 


1.5 — «EUV 


( 11 ) 


As pointed out by Becker & Bolton (20131, the local 


source approximation becomes increasingly inaccurate to¬ 
wards lower redshifts due to redshifting of Lyman limit pho¬ 


tons. Consequently, Equation (111 overestimates Fhi for a 
given Ei /,912 ■ This may be significant already at z ~ 4, since 
the quasar Hi photoionisation rate estimated from Equa¬ 


tion ( lip with our «euv = —1.7 and the rescaled emissivity 


from Masters et al. (20121 is somewhat higher than the value 
obtained by |Haardt fe Madau| ( [20l2| ) with approximately the 
same emissivity and mean free path but including redshift 
and radiative transfer effects. 

Nevertheless, Equation © provides an estimate on the 
uncertainty of Tm induced by the significant uncertainty in 
the mean EUV spectral slope (oeuv = —1.70 ± 0.61, Sec- 
tion |5.1| ). For any fixed e„,gi 2 and 2 , the ratio FHi/r H i,ref = 
3.2/ (1.5 — qeuv) quantifies the variation of the bootstrap 
realizations of Tm with respect to a reference value Fm.ref 
that uses the slope oeuv = —1.7 determined from the stack. 
As shown in Fig. |13[ the uncertain EUV slope results in a 
~ 20% uncertainty in Tm. We conclude that detailed pho¬ 
toionisation models of the IGM need to take into account 
the uncertainty in the mean EUV spectral index that is 
dominated by the large IGM correction in our study, but 
dominated by sample variance in previous work (T02, S14). 

Similarly, the uncertainty in «euv affects estimates 
of the Hen Lyman limit emissivity of quasars in models 
of Hell reionization and the post-reionization UV back¬ 
ground. Extrapolating our power-law fit to the He 11 Lyman 
limit, we obtain flux ratios /j/, 1450 / fv, 22 s = 13.8i^ 8 5 2 and 
fu, 9 i 2 / fu ,228 = lO.Olg 4 ^ 2 , respectively. Currently, the mean 
quasar SED is poorly constrained at A < 500 A due to small 


15 Note that this is an approximation due to departures of 
f(Nih,z) from a single power law (e.g. |Prochaska et al.| |2010 
Haardt &; Madau|2012). 

10 We stress that due to the low resolution we cannot fit the EUV 
continuum beneath existing EUV emission lines (S12, S14). For 
computing the photoionisation rate, this turns into an advan¬ 
tage, as we approximately account for the EUV lines that also 
contribute ionising photons. With sufficient wavelength coverage 
and spectral resolution, the SED should be integrated instead of 
the power-law continuum. 



Figure 13. Histogram of the ratio of Hi photoionisation rates 
Thi /r H i ref — 3.2/(1.5 — Qeuv)- The mean, the median, the 
16 th and the 84 th percentile are plotted as solid, dot-dashed, and 
dashed lines, respectively 


sample size (< 10 quasars) and incomplete IGM correction 
in the T02 sample. As the He II photoionisation rate depends 
on the unconstrained spectral slope at A < 228 A, models of 
the Hen-ionising background are highly uncertain. 

Lastly, we emphasize that any broken power law does 
not have any physical grounds, and thus it is a very poor 
description of the quasar SED. In principle, one should con¬ 
sider the full spectral information, although spectra do not 
usually cover energies well beyond 1 Ryd. 


6 SUMMARY AND CONCLUSIONS 

We presented the first ultraviolet average spectrum at high 
redshift ((2 em ) ~ 2.44) with proper correction for the inter¬ 
vening Lya forest and continuum absorption by employing 
the state-of the-art IGM transmission function which have 
been calibrated from multiple quasar absorption line observ¬ 
ables. The sample consists of 53 quasars selected from the 
HST survey for Lyman limit absorption systems (LLS) us¬ 
ing the Wide Field Camera 3 (WFC3) presented in O’Meara 
[etaLlpOn] ). 

The rest-frame continuum slope of the full sample shows 
a break at around 912A with a FUV spectral index ot v = 
—0.61 ± 0.01 and a softening at shorter wavelengths (a„ = 
-1.70 + 0.61). 

Our analysis highlights the fact that slope and spectral 
break in the HST composite reported by T02 are incorrect 
due to an underestimated intergalactic absorption correc¬ 
tion, especially at short wavelength where high-z quasars 
are contributing. The FUSE composite by S04 might be less 
affected by this problem given the low average redshift of 
their sample (z ~ 0.16). The COS composite by S12 (and 
S14) does not significantly differ from our spectrum within 
the uncertainties. 

Our observed broad line ratios are in good agreement 
with those predicted by the photoionisation models dis¬ 
cussed in Baskin et al. (2014), where the input quasar contin¬ 
uum has a EUV slope consistent with the one we observed. 
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The accretion disc+X-ray corona models constructed 


by Done et al. (2012) show that we need to have non zero 
values of the BH spin to match the data, with higher values 
of the BH spin as the Mbh increases (assuming that the BH 
masses derived from Civ are correct). We do not need BH 
spin if we instead consider a factor of 2 lower masses (on 
average). The latter is in agreement with previous results 
that have found BH mass values from C IV systematically 
higher once calibrated with virial mass estimators based on 
the Mg ii (e.g. Shen et al.||2008 but see also Runnoe et al. 


2013). 


Finally, we found that the lack of knowledge of the spec¬ 
tral slope results in an uncertainty on the quasar photoion¬ 
isation rate of the order of 20%. This effect should be taken 
into account by any photoionisation model. 

As a final comment, we want to stress that all previ¬ 
ous analyses were based on quasar samples selected with an 
unknown selection function (i.e. the best from the archive), 
and were characterised by small samples sizes, especially at 
short wavelengths (less than 10 AGN observations cover the 
wavelength range between 450 — 600A in the S04 sample, 
and between 10-20 AGN in the T02 sample, and only 6 con¬ 
tributing at OOOAin the S12 sample). Furthermore, errors on 
the composite were not published, so an assessment of the 
sample variance is impossible. A single statistical correction 
for the IGM absorption should not be applied for a broad 
range of redshifts and to single spectra, since LLSs cannot 
all be identified or corrected for, especially not the partial 
LLSs. 
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APPENDIX A: MULTIWAVELENGTH 
PROPERTIES OF THE WFC3 QUASAR 
SAMPLE 

We listed the multi wavelength properties of the WFC3 
quasar sample in Tables mm and|A3| 


APPENDIX B: ABSOLUTE MAGNITUDE 
RELATIONS 

Below we provide some useful relations between absolute 
magnitudes at various wavelengths. Mj(z = 2) is the mean 
absolute i-band magnitude normalised at z = 2 (Shen et 
al. 2011 ), M 1450 and M 912 are the absolute magnitudes at 
1450A and 912A respectively, estimated from our WFC3 
stack. 

Mi(z = 2 ) - Mmso = -1-28 (Bl) 


M 912 — M 1450 = 0.33 (B2) 


All these absolute magnitudes are estimated from our aver¬ 
age WFC3 spectrum rather accurately. We finally note that, 
typically, the ionising luminosity is modelled as a power law, 
thus our estimated oeuv = —1.7 can be used to evaluate the 
ionising luminosity as 


L„ = L 


912 


/ y \«EUV 

V ^912 / 


(B3) 


where the normalization, Lgi 2 , is well determined. 
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Table Al. RASS properties of the WFC3 quasar sample. 


Name 

z 

Exp. a 

s 

N H i b 

_9 

cm 

F [0.5-2]keV C 
erg s -1 cm -2 

L[0.5-2]keV d 
erg s —1 

J0755+2204 

2.319 

23963 

5.73 

4.25 ± 1.76 

1.87 ± 0.77 

J1107+6420 

2.316 

7728 

1.00 

6.74 ± 1.11 

2.95 ± 0.48 

J1119+1302 

2.394 

20105 

1.96 

9.16 ± 1.93 

4.33 ± 0.98 

J1235+6301 

2.383 

3311 

1.21 

5.56 ± 1.45 

2.61 ± 0.68 

J1253+0516 

2.398 

201 

2.17 

38.66 ± 15.21 

18.41 ± 7.24 

J1335+4542 

2.452 

616 

1.88 

10.94 ± 4.80 

5.50 ± 2.41 

J1415+3706 

2.374 

18476 

0.91 

9.47 ± 1.11 

4.40 ± 0.52 

J1454+0325 

2.368 

28496 

3.62 

24.44 ± 2.70 

11.24 ± 0.96 

J1625+2646 

2.518 

5006 

3.37 

14.40 ± 2.04 

7.73 ± 1.09 

J1651+4002 

2.343 

7649 

1.45 

17.21 ± 2.29 

7.74 ± 1.03 


[a] Exposure time. 

[b] Galactic column density (xl0 2tl ) jKalberla et al.|2005| l. 

[c] X-ray fluxes (xlO -14 ) at 0.5-2 keV calculated from the count rates in the ROSAT band using a power law spectrum with a photon 
index T = 2.0 corrected for Galactic absorption | |Kalberla et al.|2005| |. 

[d] Rest frame luminosity (xlO 45 ) in the 0.5—2 keV band (estimated from the unabsorbed X-ray flux). 


Table A2. Radio properties of the WFC3 quasar sample. 


Name 

pa 

A int 

Dist. b 

n cm 

l°g/2500 

^2500 

R f 

J0806+5041 

19.66 

0.3 

20.02 

-26.45 

-28.96 

57.08 

J0854+5327 

22.07 

0.2 

22.35 

-26.04 

-29.96 

24.52 

J1119+1302 

12.97 

0.4 

23.59 

-26.58 

-28.57 

90.48 

J1335+4542 

269.05 

0.1 

273.80 

-26.62 

-28.53 

1150.32 

J1342+6015 

1.82 

1.3 

1.83 

-26.51 

-28.76 

5.90 

J1415+3706 

409.69 

0.3 

412.16 

-26.38 

-29.05 

995.75 

J1454+0324 

4.92 

0.2 

4.95 

-26.48 

-28.79 

15.09 

J1535+4836 

107.81 

0.1 

111.14 

-26.46 

-29.04 

321.50 

J1540+4138 

18.28 

0.3 

18.77 

-26.28 

-29.45 

36.22 

J1625+2646 

9.69 

0.1 

9.95 

-26.37 

-29.08 

23.11 

J1651+4002 

43.90 

0.5 

43.96 

-26.28 

-29.28 

83.07 

J2338+1504 

45.70 

0.5 

43.75 

-26.36 

-29.15 

101.30 


[a] Integrated flux density for source in mjy. 

[b] Distance from search position to source in arcsec. 

[c] Observed radio flux density at rest-frame 6 cm in mjy. 

[d] Observed optical flux density at the rest frame 2500A in log erg s _1 cm _2 Hz _1 . 

[e] Absolute magnitude at rest-frame 2500A calculated from /2500- 

[f] Radio loudness defined as R = /6cm//2500- 
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Table A3. GALEX forced-photometry of the WFC3 quasar sample. 


Name 

z 

NUV a 

nuvJVar 

FUV a 

FUVJVar 

J0751+4245 

2.453 

0 

0 

0 

0 

J0755+2204 

2.319 

10.702 

19.899 

-0.463 

29.254 

J0806+5041 

2.457 

2.085 

0.001 

0.052 

0.679 

J0833+0815 

2.581 

0.067 

0.649 

0.291 

10.986 

J0845+0722 

2.307 

9.277 

0.155 

-0.103 

41.101 

J0850+5636 

2.464 

0.365 

19.704 

0.673 

9.939 

J0853+4456 

2.54 

0.248 

16.721 

-0.094 

27.238 

J0854+5327 

2.418 

0.734 

15.757 

-0.444 

2.087 

J0909+0415 

2.444 

-0.151 

23.423 

0.151 

9.292 

J0949+0522 

2.316 

0.351 

18.34 

-0.342 

112.288 

J1005+5705 

2.308 

0.25 

0.713 

-0.313 

61.293 

J1011+0312 

2.458 

0.252 

1.919 

0.028 

21.285 

J1053+4007 

2.482 

5.892 

0.441 

0.178 

10.585 

J1104+0246 

2.532 

2.159 

0.032 

-0.011 

2453.64 

J1107+6420 

2.316 

81.95 

4.905x 10 -5 

50.28 

3.828x 10 -5 

J1119+1302 

2.394 

0.475 

35.985 

0.044 

51.255 

J1135+4607 

2.496 

0.355 

4.004 

0.486 

12.893 

J1143+0524 

2.561 

9.468 

644.381 

0.008 

5527.879 

J1215+4248 

2.31 

3.616 

4.019 

0.015 

3909.82 

J1220+4608 

2.446 

12.665 

1.711 

14.278 

0.894 

J1228+5107 

2.45 

23.419 

1.31 

-0.002 

33.744 

J1235+6301 

2.384 

106.753 

0.107 

0.133 

35.543 

J1248+5809 

2.599 

1.221 

0.001 

-0.023 

287.568 

J1253+0516 

2.398 

-0.158 

67.674 

0.146 

43.164 

J1259+6720 

2.443 

-0.274 

4.891 

0.427 

7.997 

J1300+0556 

2.446 

0.764 

1461.025 

0.031 

5611.715 

J1302+0254 

2.414 

1.006 

10.751 

0.018 

85.795 

J1311+4531 

2.403 

3.621 

1.849 

-0.174 

7.709 

J1318+5312 

2.321 

0.372 

29.747 

0.126 

31.153 

J1323+4149 

2.44 

9.324 

0.009 

5.766 

0.022 

J1325+6634 

2.511 

0.118 

0.121 

-0.089 

499.471 

J1334+0355 

2.583 

1.219 

0.031 

-0.248 

112.13 

J1335+4542 

2.452 

1.309 

13.318 

0.459 

11.939 

J1335+4637 

2.474 

2.82 

4.541 

0.583 

6.58 

J1336+0157 

2.379 

0.378 

7.184 

0.111 

15.703 

J1342+6015 

2.399 

-0.158 

35.309 

0.077 

18.418 

J1354+5421 

2.294 

0 

0 

0 

0 

J1354+0020 

2.504 

0.134 

9.851 

-0.055 

43.665 

J1358+0505 

2.455 

0.125 

4328.555 

-0.013 

5909.892 

J1400+6430 

2.359 

7.501 

0.169 

2.517 

0.112 

J1415+3706 

2.374 

3.522 

1.408 

-0.067 

497.769 

J1454+0324 

2.368 

2.602 

11.792 

-0.168 

126.83 

J1533+3843 

2.529 

0.108 

30.193 

0.051 

22.415 

J1535+4836 

2.542 

1.847 

5.627 

0.239 

1.424 

J1540+4138 

2.516 

0 

0 

0 

0 

J1610+4423 

2.588 

1.151 

6.765 

-0.243 

13.434 

J1625+2943 

2.357 

-0.054 

12.542 

0.233 

1.771 

J1625+2646 

2.518 

57.894 

0.382 

0.267 

15.956 

J1651+4002 

2.343 

0.184 

11.125 

0.189 

108.886 

J1724+5314 

2.547 

0.608 

63.628 

-0.206 

230.608 

J2111+0024 

2.325 

0 

0 

0 

0 

J2136+1029 

2.555 

0.276 

0.643 

0.137 

6.969 

J2338+1504 

2.419 

4.765 

0.101 

0.435 

0.234 


[a] Fluxes in AB nanomaggies (monochromatic AB fluxes are defined as Fab = Fab, nano X 10 11 tx0.4 

[b] Pipeline inverse variance (IVAR) in AB nanomaggies. 
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